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We propose and develop an efficient implementation of the robust tabu search heuristic for sparse quadratic 
assignment problems. The traditional implementation of the heuristic applicable to all quadratic assignment 
problems is of 0{N 2 ) complexity per iteration for problems of size N . Using multiple priority queues to 
determine the next best move instead of scanning all possible moves, and using adjacency lists to minimize 
the operations needed to determine the cost of moves, we reduce the asymptotic (N — >• oo) complexity per 
iteration to O(NlogN). For practical sized problems, the complexity is O(N). 
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1. Introduction 

The quadratic assignment problem (QAP) is 
a combinatorial optimization problem first intro- 
duced by Koopmans and Beckman (1957). It is 
NP-hard and is considered to be one of the most 
difficult problems to be be solved optimally. The 
problem was defined in the following context: A 
set of N facilities are to be located at N locations. 
The distance between locations i and j is Dij and 
the quantity of materials which flow between facili- 
ties i and j is Fij. The problem is to assign to each 
location a single facility so as to minimize the cost 
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where p(i) represents the location to which facility 
% is assigned. It will be helpful to think of the N 
facilities and the matrix of flows between them in 
graph theoretic terms as a graph of N nodes and 
weighted edges, respectively. 

There is an extensive literature which addresses 
the QAP and is reviewed in Pardalos et al. (1994), 
Cela (1998), Anstreicher (2003), Loiola et al. (2007, 
and James et al. (2009a). With the exception 
of specially constructed cases, optimal algorithms 
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have solved only relatively small instances with 
N < 36. Various heuristic approaches have been 
developed and applied to problems typically of size 
N « 100 or less. One of the most successful heuris- 
tics to date for large instances is robust tabu search, 
RTS, (Taillard (1991). The use of tabu search for 
the quadratic assignment problem has been stud- 
ied extensively (Drezner (2005), Hasegawa et al. 
(2000), James et al, (2009a, 2009b), McLoughlin 
and Cedeno (2005), Misevicius (2007), Misevicius 
and Ostreika (2007), Skorinkapov (1994), and Wang 
(2007)). Some of the best available algorithms for 
the solution of the QAP are the hybrid genetic al- 
gorithms that use tabu search as an improvement 
mechanism. (See Drezner (2008)). 

Here we will consider the robust tabu heuristic 
applied to sparse QAP instances. That is, the num- 
ber of non-zero entries in the either the flow matrix 
and/or the distance matrix is of O(N) as opposed 
to 0(N 2 ). Without loss of generality we will as- 
sume the flow matrix is sparse. Many real world 
problems are sparse. In fact, this work was moti- 
vated by the study of random regular sparse graphs. 
These graphs are very robust to partitioning and 
collapse due to removal of nodes or edges. We are 
interested in the problem of determining how to as- 
sign the nodes of such a graph to locations in a 
metric space such that the total edge length of the 
graph is minimized; this problem maps directly to 
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a quadratic assignment problem. 

There has been some previous work on sparse 
quadratic assignment problems. Milos and Magirou 
(1995) developed a Lagrangian-relaxation lower- 
bound algorithm for sparse problems and Panos 
et al. (1997) developed a version of their GRASP 
heuristic for sparse problems. However, to the best 
of our knowledge, an efficient implementation of the 
robust tabu heuristic for sparse QAP instances has 
not been proposed. 

2. Background - the tabu heuristic 

The tabu heuristic for the quadratic assignment 
problem consists of repeatedly swapping locations 
of two nodes. A single iteration of the heuristic 
consists of 

(a) Determining the move which most decreases 
the total cost. Under certain conditions (see 
Section 4), if a move which lowers the cost is 
not available, a move which raises the cost is 
made. So that cycles of the same moves are 
avoided, the same move is forbidden (taboo) 
until a specified later iteration; we call this 
later iteration the eligible iteration for a given 
move. This eligible iteration is traditionally 
stored in a tabu list or tabu table. 

(b) Making the move. 

(c) Recalculating the new cost of all moves. 

The process is repeated for a specified number 
of iterations. Traditional implementations of ro- 
bust tabu search require 0(N 2 ) operations per it- 
eration. The complexity of 0(N 2 ) is achieved by 
maintaining a matrix containing the cost A(p, u, v) 
of swapping u and v for all u and v, given a current 
assignment p. 

The complexity of the each step above is as fol- 
lows: 

(a) 0(N 2 ) - all possible N(N - l)/2 moves are 
considered. The cost of each move is retrieved 
from A(p, r, s) 

(b) O(l) - the locations of the two swapped nodes 
are simply transposed. 

(c) 0(N 2 ) - based on the following observations of 
Taillard (1991): 



Following Taillard (1991), starting from an as- 
signment of facilities p let the resulting assign- 
ment after swapping facilities r and s be p'. 
That is: 

p (k) — p(k) k^r,s 
p'(r) = p(s) 

p'(s) = p(r). (2) 

For a symmetrical matrix with a null-diagonal, 
the cost A(p 1 r, s) of swapping r and s is: 

N N 

A(p,r,s) = ^^(f,jfl p (,),p( 3 ) ~ F i,3 D p'(i), P 'U)) 

i=l j=l 

= 2 ( Fs ' k ~ ^fe)( £ 'p(s),P(fc)--Dp(r),P(^) 

k^r.s 

To calculate A(p',u, v) for any u and v with 
complexity O(N) , we can use equation xx. 
For asymmetric matices or matrices with non- 
null diagonals, a slightly more complicated ver- 
sion of Equation ^ also of complexity O(N) 
is given by Burkhard and Rendl (1984) . 

To calculate A(p',u, v) in the case that the 
swapped facilities u and v are different from 
r or s, we use the value A(p, u, v) calculated in 
the previous iteration and find: 

A(p', u, v) = A{p, u, v) + 2(F r> „ - F r , v + F s , v - F s>u ) 

' {J-^p 1 (s) 5 p' (u) ^p'(s),p'(v) -^p' (r) ,p' (v) Dpi ( r ) pi (w)) (^ 

(i) the cost of moves which do not involve 
the two nodes in the previous move can be 
calculated in time O(l). There are 0(N 2 ) 
of these moves. 

(ii) The cost of moves which do involve the 
two nodes in the previous move must be 
calculated from scratch. There are 0(N) 
of these moves and the complexity of cal- 
culating each is O(N). 

3. Approach 

To reduce the complexity of step (a), instead of 
scanning all possible moves, we use multiple prior- 
ity queues (PQs) to determine the best move. A 
priority queue is a data structure for maintaining 
a set of elements each of which has an associated 
value (priority). A PQ supports the following op- 
erations: 
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• Insert an item 

• Remove an item 

• Return the item with the highest value 

Priority queues are used to efficiently find an item 
with the highest value without searching through 
all of the items. 

The maximum complexity of PQ operations is 
0{logN). We will see below that there will be O(N) 
insertions and deletions in the PQs for each itera- 
tion so the asymptotic complexity of this step is 
reduced to O(NlogN). Furthermore, we will show 
that for problems of any practical size, PQ opera- 
tions are not the determinant of total complexity. 

The complexity to recalculate the cost of moves 
in step (c), can be reduced to 0{N) as follows: 

• As in the traditional robust tabu implementa- 
tion, the cost of moves which do not involve the 
two nodes in the previous move can be calcu- 
lated in time O(l). On average, there are 2(k) 
nodes which are connected to the two nodes in 
the previous moves, where (k) is the average 
degree (average number of nodes adjacent to a 
given node) of the graph corresponding to the 
flow matrix. For each of these 2{k) nodes we 
must calculate the cost of N— 1 possible moves. 
Thus, the cost is 0((k)N). 

• The cost of moves which do involve the two 
nodes in the previous move must be calcu- 
lated from scratch. There are O(N) of these 
moves and the complexity of calculating each 
is 0((k)) since the cost of a node, u, being 
in a specific location depends only on the on- 
average k nodes adjacent to u. 

Thus the complexity of step (c) is reduced to O(N). 



4. Implementation 

To describe our implementation, we must first 
describe the rules for determining the next move of 
Taillard's robust tabu heuristic (Taillard (1991)). 
The following definitions for the possible state of a 
potential move are useful: 

(i) If the current iteration is less than or equal to 
the eligible iteration, the move is ineligible. 

(ii) If the current iteration is greater than the eli- 
gible iteration, the move is authorized. 




Figure 1: Priority queues used in the implementation. 

(iii) If the current iteration minus an aspiration 
constant is greater than the eligible iteration 
the move is aspired. 

The rules for determining the next move can then 
be stated as (Taillard (1991)): 

(1) If a move which decreases the lowest total cost 
found so far is available, the move which most 
decreases this total cost is chosen, independent 
of whether the move is ineligible, authorized or 
aspired. 

(2) If no move meets criterion (1), the aspired 
move, if one is available, which most decreases 
the current total cost is chosen. 

(3) If no moves meet criteria (1) or (2), the lowest 
cost authorized move is chosen. 

To implement these rules for sparse problems, we 
use two types of PQs: delta PQs which contain the 
cost delta for a given move and tabu PQs which 
contain entries ordered by the eligible iteration for 
the move. The tabu PQs control the change of state 
of a move. The delta PQs determine the lowest cost 
move in each state. Five PQs are used: 

• ineligible tabu PQ - This PQ contains moves, 
ordered by eligible iteration, which are in the 
ineligible state. This PQ allows us to efficiently 
determine when the state of a move can be 
changed to authorized. 

• authorized tabu PQ - This PQ contains moves, 
ordered by eligible iteration, which are in the 
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authorized state. This PQ allows us to effi- 
ciently determine when the state of a move can 
be changed to aspired. 

• ineligible delta PQ - This PQ contains moves, 
ordered by the cost of the move, which are in 
the ineligible state. This PQ together with the 
two other delta PQs allows for efficient deter- 
mination of the overall lowest cost move as re- 
quired by rule 1. 

• aspired delta PQ - This PQ contains moves, 
ordered by the cost of the move, which are in 
the aspired state. This PQ allows for efficient 
determination of the lowest cost aspired move 
as required by rule 2. 

• authorized delta PQ - This PQ contains moves, 
ordered by the cost of the move, which are in 
the authorized state. This PQ allows us to 
determine the lowest cost authorized move as 
needed by rule 3. 

As illustrated in Fig. [TJ moves are inserted and 
removed in the PQs under the following circum- 
stances: 

• At initialization all moves are inserted into the 
ineligible PQs. 

• At the beginning of each iteration, any moves 
on the ineligible PQs which become authorized, 
because the iteration has increased by one, are 
moved from the ineligible PQs to the autho- 
rized PQs. 

• At the beginning of each iteration, any moves 
on the authorized PQs which become aspired, 
because the iteration has increased by one, are 
moved to the aspired PQ. 

• After each move, any move for which the move 
cost or eligible iteration has changed is re- 
moved from the PQ in which it is present and 
inserted in the appropriate PQ based on move 
state and move cost. (However, see lazy up- 
date discussion below.) 

Using these PQs we obtain exactly the same results 
as the traditional robust tabu implementation. 

5. Lazy Update 

We minimize the time updating PQs by perform- 
ing lazy updates. After a change in the eligible 



iteration or move cost, if the state is changed or 
the value is increased, we update the PQs involved; 
otherwise, we perform a lazy update and store the 
value in a data structure associated with the move 
and only do the update in the PQ when and if the 
move becomes the move with the smallest value in 
the PQ. This use of lazy updates significantly de- 
creases the time spent on PQ operations. 




1 2 1 3 

N 

Figure 2: The ratio of the time per iteration for the non- 
sparse implementation to the time per iteration for the sparse 
implementation versus the problem size, N. The slope of the 
plot is approximately 1.0 which is consistent with a reduction 
of O(N) of the computational complexity 




N 

Figure 3: The time per iteration for the non-sparse (upper 
plot) and sparse (lower plot) implementations. 



6. Numerical Results 

We test our algorithm on instances with N lo- 
cations on a square grid with a Euclidean metric. 
The flow matrix for the N facilities corresponds to 



4 




k 

Figure 4: The time per iteration using the sparse implemen- 
tation for N = 400 versus degree k. 

an adjacency matrix for a k- regular random graph; 
that is, each facility has flows of value one to k other 
random facilities. We run both the traditional non- 
sparse implementation and our new sparse imple- 
mentation. The non-sparse implementation is the 
code of Taillard (1991). To implement the priority 
queues in the sparse implementation, we use the 
complete balanced tree code of Marin and Cordero 
(1995) . We perform a minimum of 10 4 iterations 
for each instance. 

For N — 2500, priority queue update operations 
consume just 2% of the total time. Assuming a 
logN behavior for PQ operations, it is not until N 
becomes astronomically large that the PQ opera- 
tions would take the majority of the time and the 
behavior of our implementation crosses over from 
0(N) to O(NlogN) complexity. We thus expect 
O(N) complexity for our sparse implementation for 
any practical sized problems. 

In Fig. [3J for k = 3, we plot the ratio r of the 
time per iteration of the non-sparse implementation 
to the time per iteration of the sparse implementa- 
tion versus N. As expected, the plot is consistent 
with r ~ N x where x ss 1.0. This reflects the 
reduction in complexity from 0(N 2 ) to O(N). In 
Fig. |3] we plot separately the time per iteration for 
the original and the sparse implementations. Con- 
sistent with Fig. the slopes on the log-log plot 
differ by 1 but the slopes are 2.5 and 1.5, respec- 
tively, as opposed to the theoretical values of 2.0 
and 1.0. As explained in Paul (2007) and Saavedra 
and Smith (1995), this is due to the finite size of 
processor cache memory; as the problem size (and 
memory needed) increases, there are a smaller per- 
centage of cache hits causing slower operation. 



In Fig. 3] we plot the time per iteration versus 
values of degree k. The plot is linear with a slight 
deviation with increasing k. This deviation from 
linear is due to the fact that as k increases there 
is an increasing chance that a node u will be con- 
nected to both nodes involved in the previous move; 
however, the updated costs of moving node u must 
be calculated only once. 

We also performed numerical experiments on 
a class of problems known in the literature (see 
Drezner et al. (2005)). These problems denoted as 
DRExx are sparse and are designed to be very dif- 
ficult to solve with heuristics. We obtained results 
similar to those for described above; performance is 
O(N). 

7. Summary 

For sparse quadratic assignment problems, we re- 
duce the asymptotic complexity per iteration of ro- 
bust tabu search to O(NlogN) from 0(N 2 ); im- 
practical size problems, the complexity is reduced 
to O(N). Central to achieving this reduction is the 
use of multiple priority queues and lazy updates to 
these queues. The code which implements our ap- 
proach and test QAP instances used for this paper 
are available as supplementary material. 
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